﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace RouteConvertTool
{
    public static class UtmConverter
    {
        public static void ConvertLatLonToUtm(double latitude, double longitude, out double easting, out double northing, out int zoneNumber)
        {
            // WGS84椭球参数
            const double a = 6378137.0;              // 长半轴
            const double f = 1.0 / 298.257223563;   // 扁率
            const double k0 = 0.9996;               // 比例因子
            const double falseEasting = 500000.0;
            const double falseNorthing = 10000000.0;

            // 计算UTM带号 (1-60)
            zoneNumber = (int)((longitude + 180) / 6) + 1;

            // 计算中央经线 (单位：度)
            double centralMeridian = (zoneNumber - 1) * 6 - 180 + 3;

            // 转换为弧度
            double latRad = ToRadians(latitude);
            double lonRad = ToRadians(longitude);
            double centralMeridianRad = ToRadians(centralMeridian);

            // 椭球参数计算
            double e = Math.Sqrt(2 * f - f * f);     // 偏心率
            double e2 = e * e / (1 - e * e);         // 第二偏心率平方
            double n = a / Math.Sqrt(1 - Math.Pow(e * Math.Sin(latRad), 2));
            double T = Math.Pow(Math.Tan(latRad), 2);
            double C = e2 * Math.Pow(Math.Cos(latRad), 2);
            double A = (lonRad - centralMeridianRad) * Math.Cos(latRad);

            // 子午线弧长计算
            double M = a * ((1 - e * e / 4 - 3 * Math.Pow(e, 4) / 64 - 5 * Math.Pow(e, 6) / 256) * latRad
                       - (3 * e * e / 8 + 3 * Math.Pow(e, 4) / 32 + 45 * Math.Pow(e, 6) / 1024) * Math.Sin(2 * latRad)
                       + (15 * Math.Pow(e, 4) / 256 + 45 * Math.Pow(e, 6) / 1024) * Math.Sin(4 * latRad)
                       - (35 * Math.Pow(e, 6) / 3072) * Math.Sin(6 * latRad));

            // 计算东偏移(Easting)
            double eastingValue = k0 * n * (A + (1 - T + C) * Math.Pow(A, 3) / 6
                                  + (5 - 18 * T + T * T + 72 * C - 58 * e2) * Math.Pow(A, 5) / 120)
                                  + falseEasting;

            // 计算北偏移(Northing)
            double northingValue = k0 * (M + n * Math.Tan(latRad) * (A * A / 2
                                      + (5 - T + 9 * C + 4 * C * C) * Math.Pow(A, 4) / 24
                                      + (61 - 58 * T + T * T + 600 * C - 330 * e2) * Math.Pow(A, 6) / 720));

            // 南半球处理
            if (latitude < 0)
                northingValue += falseNorthing;

            easting = Math.Round(eastingValue, 6);
            northing = Math.Round(northingValue, 6);
        }

        private static double ToRadians(double degrees)
        {
            return degrees * Math.PI / 180.0;
        }
    }

}
